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Abstract 

We present a hybrid method for the simulation of colloidal systems, that combines molecular 
dynamics (MD) with the Lattice-Boltzmann (LB) scheme. The LB method is used as a model 
for the solvent in order to take into account the hydrodynamic mass and momentum transport 
through the solvent. The colloidal particles are propagated via MD and they are coupled to the 
LB fluid by viscous forces. With respect to the LB fluid, the colloids are represented by uniformly 
distributed points on a sphere. Each such point (with a velocity V(r) at any off-lattice position 
r) is interacting with the neighboring eight LB nodes by a frictional force F = £o(V(r) — u(r)) 
with £o being a friction force and u(r) being the velocity of the fluid at the position r. Thermal 
fluctuations are introduced in the framework of fluctuating hydrodynamics. This coupling scheme 
has been proposed recently for polymer systems by Ahlrichs and Diinweg [J. Chem. Phys. Ill, 
8225 (1999)]. We investigate several properties of a single colloidal particle in a LB fluid, namely the 
effective Stokes friction and long time tails in the autocorrelation functions for the translational and 
rotational velocity. Moreover, a charged colloidal system is considered consisting of a macroion, 
counterions and coions that are coupled to a LB fluid. We study the behavior of the ions in a 
constant electric field. In particular, an estimate of the effective charge of the macroion is yielded 
from the number of counterions that move with the macroion in the direction of the electric field. 

PACS numbers: 47.65. +a, 82.70.Dd, 47.11. +j, 83.85.Pt, 05.20.Dd, 07.05.Tp, 66.20.+d 
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I. INTRODUCTION 



In colloidal suspensions, big molecules (colloids) with a typical size of the order of 10 nm 
to several /im are immersed in an atomistic solventA* 2 -. A detailed modeling of such systems 
would be provided by a molecular dynamics (MD) simulation where both colloids and solvent 
atoms are propagated via Newton's equations of motion. However, this approach has severe 
drawbacks due to the large size difference between colloids and solvent particles. One would 
have to take into account the microscopic details of a large amount of solvent particles which 
are irrelevant on the typical length and time scale of the colloids. In order to circumvent this 
problem, one may describe the interactions between the colloids by an effective potential, 
thus avoiding the explicit consideration of the solvent's degrees of freedom^. But then the 
hydrodynamic mass and momentum transport through the solvent is completely neglected 
and hence, hydrodynamic interactions between colloidal particles are not taken into account. 
But it is well-known that many transport properties in colloidal suspensions are affected by 
hydrodynamic interactions^. 

In recent years many efforts have been undertaken to model colloidal suspensions by 
mesoscopic simulation techniques. The idea is to describe the solvent on a coarse-grained 
level, whereby the properties of the solvent on hydrodynamic time and length scales are 
correctly recovered. Different approaches exist in the literature that range from particle- 
based methods such as dissipative particle dynamics (DPD) 3 ' 4 ' 5 ' 6 and stochastic rotation 
dynamics (SR-D) 7 ! 8 ! 9 ' 10 ' 11 to the Lattice-Boltzmann (LB) metho d 12113114 where the Navier- 
Stokes equations are solved on a lattice via a kinetic equation. All these methods have their 
advantages and disadvantages, and they might provide complementary information for a 
given problem. 

In this paper, we present a hybrid MD/LB method for the simulation of colloidal systems. 
In this approach, the coupling between colloids and solvent (LB fluid) is realized locally at 
the surface of the colloids. The modeling of the solvent by a LB fluid has the advantage that 
hydrodynamic properties of the solvent (e.g., the shear viscosity rj) are incorporated as input 
parameters, and thus they can be varied over a broad range (compared to particle-based 
methods). Moreover, the LB method allows to implement easily the rotational degrees of 
freedom of the colloids, an issue that is difficult to handle in particle-based models of the 
solvent. 
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Several methods have been proposed for the introduction of colloidal particles in a LB 
fluid. In Ladd's metho d 12 ' 13 ' 15 , a colloid is represented by its surface which cuts some of 
the links between lattice nodes of the LB fluid. Halfway along these links, boundary points 
are placed where the one-particle distribution functions, that represent fluid populations, 
are bounced back such that no-slip boundary conditions are obtained. A different approach 
for the particle-fluid coupling has been proposed by Lobaskin et alM Here, a colloid is 
represented by point particles that are connected with each other by springs (modeled by 
a FENE potential). The latter point particles interact with the fluid locally via a friction 
f orcQ i7 1 i8 1 i2 COU p]j n g scheme that we present in this paper can be seen as one in 

between Ladd's bounce back rules and the frictional coupling scheme of Lobaskin et al. As 
in Ladd's method, each colloid is represented by boundary points on its surface and the 
total force and torque that is exerted on a colloid by the fluid is obtained by summing up 
all the contributions from the boundary points. Different from Ladd's method we use local 
friction forces between the boundary points and the fluid as proposed by Lobaskin et al. 
However, no interaction potential between boundary points is required in our scheme. We 
also avoid some inherent problems of Ladd's method. Our scheme does not introduce shape 
fluctuations of a moving colloidal particle as well as mass fluctuations in the fluid that are 
induced by the bounce back rules on a moving colloidal particle. 

Many attempts have been devoted very recently to the development of simulation methods 
for charged colloids, using the LB method or finite difference schemes to model hydrody- 
namic interactions. These approaches consider charged colloids in the framework of the 
primitive model-, i.e. as a system of negatively charged macroions, small counterions of pos- 
itive charge and small coions of negative charge. In Refs .^21*22*2^ the small ions are modeled 
as charge densities on the level of the Poisson-Boltzmann equation. The propagation of 
these charge densities is achieved by the so-called electrokinetic equations^ that couple the 
Navier-Stokes equations with a dynamic generalization of the Poisson-Boltzmann equation. 
The description of charged colloids by the electrokinetic equations has some drawbacks: Of 
course, correlations between the ions are neglected, since the Poisson-Boltzmann equation 
is a mean-field description. Moreover, it is not clear how one can properly introduce ther- 
mal fluctuation in the electrokinetic equations and thus, so far only calculations at zero 
temperature are possible. Therefore, Lobaskin et alM have followed a different strategy by 
considering explicitly small ions that are moved together with macroions by their hybrid 
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MD/LB scheme. Using MD/LB approach, we also model charged colloidal systems in this 
way, and, as an application we consider below the problem of a charged colloidal system in 
an external electric field. 

The rest of the paper is organized as follows: In the next section, we give a brief intro- 
duction to the LB method. Then, we present our scheme to couple colloidal particles with 
a LB fluid (Sec. Ill) and the use of this scheme in a hybrid MD/LB method (Sec. IV). 
Applications are shown in Sec. V where a neutral colloid in a LB fluid is considered as well 
as a macroion and small ions in a LB fluid that are subject to an external electric field. 
Finally, we summarize and discuss the method and the results. 

II. LATTICE BOLTZMANN METHOD 

In this section, we briefly describe the LB method used in our work. A detailed description 
can be found in several review article a 12 ' 1 ? and in the book by Succi^. 

In the LB method a discretized version of a simple kinetic equation is solved numerically 
on a lattice. The central quantity in this equation is an one-particle distribution function 
rii(r,t) which gives the particle population on a lattice node r at time t with a discrete 
velocity Cj. The discrete space of velocities {cj} has to be constructed such that the un- 
derlying kinetic equation is consistent with the Navier-Stokes equations. In particular, the 
{cj} should not introduce any artificial anisotropies in the hydrodynamic equations. Several 
choices for {cj} are possible. In this work, the velocity space consists of 18 vectors placed on 
lattice nodes of a cubic lattice, 6 of these point to the nearest neighboring nodes and 12 to 
the next-nearest neighboring nodes. These velocity vectors have the absolute values 1 a/r 
and v2 a/r, respectively, with a being the lattice spacing and r the elementary time unit. 
Note that the latter velocity space can be constructed by projecting the 24 unit vectors of 
a four dimensional FCHC lattice onto three dimensions. 

The equations of motion for the rii(r,t) can be written as 

riiir + dT,t + t) = rii(r,t) + Ai[{m(r, *)}], (1) 

where the collision operator A, describes the change of rij due to collisions between the 
particles (of course Aj has to be further specified). The moments of rij(r, t) are hydrodynamic 
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fields, namely the mass density 

18 



p(r,t) = ^TK(r,t) , (2) 



the momentum density 



and the momentum flux 



i=l 
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j( r ' t ) = ^2 n i( r >t) c i > ( 3 ) 



i=l 
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II(r,t) = ^ni(r,t)ciCi . (4) 

i=l 

Note that j can be also written as j = pu with u denoting the flow velocity. 

The aim is to construct the collision operator Aj in as simple a form as possible such 
that the time evolution for the moments in Eqs. (fl|)— (|4*|) is consistent with the Navier-Stokes 
equations. To this end, we assume that Aj describes small deviations from equilibrium. 
Hence, it can be given in a linearized form, 

18 

A,[K(r,t)}] = ^L,,K-(r,t) - n^) , (5) 
i=i 

where n eq is the local equilibrium distribution function and denotes a matrix element of 
the collision operator. 

The equilibrium population n eq can be expressed a o 12 > 14 

1 1 



n e - q 



p + -^pu ■ c, + — puu : (CjCj - c 2 s l) 



(6) 



where n eq depends only on p and u. The parameter c s in Eq. © denotes the sound velocity. 
The values for c s as well as for the coefficients a Ci can be fixed by a Chapman-Enskog 
expansion such that the Navier-Stokes equations in the limit of small u are recovered. This 
yields c s = a/1/2. The parameters a Cl can have the two values a 1 and where the 
superscript corresponds to the absolute value of the velocity vector. In our case, we have 
a 1 = 1/12 and = 1/24. The conservation of mass and momentum requires that p 
and j are unaffected by the collision operator. The kinematic viscosity v is related to the 
eigenvalue A of that corresponds to the eigenvector c^c^ (herein, the indices a, (3 with 
a 7^ (3 denote the Cartesian components): 



In the following we choose A = —1.75 which corresponds to v = 0.0238 a 2 /r. 

All non-hydro dynamic modes (ghost modes), that are due to the use of a lattice, are 
suppressed by setting the corresponding eigenvalues of to -1. Since we are interested in 
incompressible LB fluids also the eigenvalue related to the bulk viscosity is chosen to be -1. 

As shown by Ladd^, the LB scheme as defined by Eqs. (fljl— (|7J) allows also the introduction 
of thermal fluctuations. This can be done in the framework of fluctuating hydrodynamics^. 
To this end, one adds an additional stochastic term in Eq. (0) which represents the 
addition of thermal fluctuations to the stress tensor, 

n 'i = --jzL^^f 3 ■ ( 8 ) 

s a/3 

The random stresses a' a/3 are random numbers with zero mean and white noise behavior, i.e. 

= A5 rr >5tt> (fioryfipS + £a<5<5/3 7 — -5q/3^ 7 5^ (9) 

where the 5's are Kronecker deltas and A is a constant which has to be chosen such that 
the fluctuation-dissipation theorem holds. The latter condition requires the choice A = 
2rjk B T\ 2 . 

An elementary step in a LB simulation can be decomposed in a collision step and a 
propagation step. In the collision step, the interaction between the "particles" at node r at 
the collision time t* is taken into account that results in the postcollision function 

<(r, **) = m(r, t*) + At[{m}] + r4(r, **)■ (10) 

Using the n*(r,t*), the n« can be updated in the propagation step, 

n i (r + c i ,t + l) = <(r,t*) . (11) 

In our implementation we have omitted non-linear terms in u and thus, we describe a fluid 
on the level of the linearized Navier-Stokes equations (for more details see Refs.— *^) . 

III. THE COUPLING OF THE LB FLUID TO COLLOIDAL PARTICLES 

Consider a Brownian particle of mass M in a solvent. In first approximation, one can 
describe the motion of the particle without taking into account explicitly the solvent par- 
ticles. To this end, one assumes that, on the typical Brownian time scale, the collisions of 
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the particle with the solvent particles can be modeled by Gaussian random forces f r . These 
forces lead to a systematic friction force — £oV(t) on the particle where £o is the friction 
coefficient and V(£) the velocity of the particle at time t. The resulting equation of motion 
is a Langevin equation (see, e.g., Ref.— ), 

M^ = F c -£ V(t) + f r . (12) 

Here, r is the position of the particle and F c denotes a conservative force due to the interac- 
tion with other particles. The Cartesian components of the forces, f^ a , are random numbers 
that are uncorrelated with zero mean, i.e. 

(f t , a (r,t)) = (13) 
(frATMA^J)) = AS ^S(r - r')5(t - . (14) 

The amplitude A is given by the fluctuation-dissipation theorem, A = 2A;bT£ . 

In Eq. ()12|) the interaction of the Brownian particle with the solvent is described by the 
force Fb = — £oV(t) + f r . By this force, transport of momentum in the fluid due to the 
motion of the particle is not taken into account and thus one does not recover the correct 
hydrodynamic behavior of the Brownian particle. This is reflected, e.g., in the behavior 
of the velocity autocorrelation function C v (t). For a single particle, propagated according 
to Eq. fT2"|) (with F c = 0), C v (t) decays exponentially, C v (t) oc exp(— £ t/M), whereas the 
correct hydrodynamic behavior is a power law decay at long times, C v (t) oc t~ 3 / 2 (the so- 
called long-time tail)^. However, in order to incorporate hydrodynamics, Eq. (|T2*j) can be 
easily modified by coupling the Brownian particle to a LB fluid. The essential point is to 
replace the absolute velocity V(t) of the particle in the frictional force term by its velocity 
relative to the fluid. The force Fb is then modified to 

F B = -MVM)-u(r,t)) + f r (15) 

where u(r, t) is the velocity of the fluid at the position of the particle. Since the position 
of the particle is continuous in space, we have to use an interpolation scheme to determine 
u(r,t). As in Ref . 18 ' 19 , we use a linear interpolation scheme that estimates u(r,t) from the 
eight nearest lattice nodes around r. Of course, Newton's third law requires that whenever 
the Brownian particle is subject to a force F B as given by Eq. (fTo]). the force — F B has to 
be applied to the fluid nodes with which the particle interacts (see Refs . 18 ' 19 ). One may 
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wonder why fluctuations have to be added both to the LB fluid and the Brownian particle. 
However, it was shown in Ref.~ that only then the fluctuation-dissipation theorem holds 
for the total system of a particle in a LB fluid coupled via Eq. (|15|). 

Up to now we have considered the Brownian particle as a point particle which has, in 
particular, no rotational degrees of freedom (however, its effective hydrodynamic radius 
with respect to the LB fluid is of the order of the lattice spacing a). In order to model an 
extended object, such as a sphere of radius R, the force coupling as given by Eq. ([To]) can 
be generalized as follows: As in Ladd's metho d 12 i 13 the sphere is represented by boundary 
points on its surface, whereby the surface is permeable for the LB fluid. In our method, the 
boundary points are placed uniformly on the surface of the sphere. The uniform distribution 
of boundary points can be achieved by the following iterations: One starts from points that 
are placed at the corners of an octahedron. The envelope of this octahedron is a sphere 
of radius R, the centre of which is denoted by R cm in the following. Then one places new 
points halfway along the twelve edges of the octahedron and shifts these points such that 
they sit on the surface of the latter sphere. As a result, a polygon with 18 corners and 
48 edges is obtained. In this object, again new points are placed halfway along the edges 
and shifted to the surface of the enveloping sphere. This object has now 66 corners, i.e. 66 
boundary points, and is used in this work as a model for a colloidal particle. A snapshot is 
shown in Fig. We consider that each of the boundary nodes has a mass M/66 where M 
is the mass of the colloidal particle. Each of the boundary points is coupled to the LB fluid 
according to Eq. (fT5|) . The total force on the particle is then determined by the sum 

66 

F B ,tot = ^F B (r lb ) (16) 

i b =l 

where r ih denotes the positions of the 66 boundary points. The torque that is exerted on 
the particle by the fluid is updated via 

66 

T B = ^F B (r,, b )x(r ib -R cm ). (17) 

i b =l 

F B ,tot and T B can be used to update the translational and the rotational velocity of the 
colloidal particle which we denote by V and ujq, respectively. The simplest discretized form 
of the equations of motion is the Euler algorithm which is of first order in the time step h 
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for the integration, 

Vo(t + h) = V (t)+F B>tot — (18) 
u{t + h) = uj(t) + T B j . (19) 

In Eq. f!19|) . / is the moment of inertia of the particle. The Euler algorithm has of course 
very bad properties with respect to stability, temperature drift etc. In the next section, we 
present an integrator based on the Heun algorithm 28 that is of second order in h. 

The boundary points that we use for the representation of the surface of the spherical 
particle are not rotated according to the updated rotational velocity uj . These points are 
fixed with respect to the centre of mass position R cm of the particle. The idea is that 
the surface of a spherical particle can be always represented by the same set of uniformly 
distributed points if the surface grid that they form is fine enough. We have checked that 
the sphere with 66 boundary points (see Fig.P) is sufficient for particle radii up to R = 5.0a. 
In this case, no change in any physical properties is seen if the particle is decorated with 
more boundary points. 



IV. THE HYBRID MD/LB SCHEME 

Now, we discuss the algorithm that we use for the integration of the equations of motion 
of a system of colloidal particles coupled to a LB fluid. To this end, we use a generalized 
velocity Verlet algorithm 29 . For the colloidal particles we use the model with 66 boundary 
points that we have described in the previous section. 

The force on particle i can be decomposed into three terms, 

F t ot,i = F Ci j + Ff 5 j + F ri j, (20) 

where F C) j is the contribution from conservative forces. Ff j and F r) j denote friction and 
random forces, respectively. To compute Ff^, one sums up the contributions from the 66 
boundary points, 

66 

Ff,i = ~J2 fob t v i + U X fob - r i) - u ( r J] > ( 21 ) 

»b=l 

where and Vj are respectively the position and the velocity of the centre of mass of particle 
i, r ih is the position of the boundary node of particle i, uj is its angular velocity, and the 
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fluid velocities at the boundary points are denoted by u(r ib ). £ ,j b is the friction coefficient 
of boundary point % and thus the total friction coefficient £o is given by £o = J2i b £o,i b • In 
the following, each of the 66 £ ,i b is assigned the same value £0/66 for a given particle. 
The random force can be written as 

66 

F r ,« = v / 2vrA; B T ^ 6 ih . (22) 

»b=i 

Here, the 6i is a vector of random numbers with the Cartesian components 0j jQ (a G {x, y, z}) 
for which 

(M*)) = ( 23 ) 
(fi, a (i)^(0) = 5«M(*-0 (24) 

holds. Note that it is not necessary to use random numbers #j )Qi with a Gaussian distribution 
for the numerical integration. It was shown in Ref.— that it is sufficient to use uniform 
random numbers that fulfill the requirements as given by Eqs. f!2^l) and ()24j) . 

The update of the centre of mass positions of the particles is similar to that in the velocity 
Verlet algorithm 29 . In an integration time step h the position change from i"j(0) to 

Ti(h) = r«(0) + hvi(0) 

+^(F c , i (0) + F M (0)) + F r , i (/ i ). (25) 

For the determination of the velocities Vi(h), the friction force Ff j at time t = h is needed 
if we want to apply the velocity Verlet scheme. But in this case, the problem arises that 
Ff,j(/i) depends itself on Vi(h). This problem can be solved if we first approximate Vi(h) in 
an Euler step to obtain 

v*(h) = Vi (0) + A (F Cji (0) + F M (0)) + ^F Tji (h) , (26) 

where the star indicates that we use this velocity only for an estimate of Ff^/i). Also the 
angular velocity u>i(h) is updated by an Euler step, 

h 66 

Wi{h) =Ui(0) + jJ2 F f (0) x [r, b - n(h)] . (27) 
»b=i 

With Eqs. and (|2*7j) we yield the friction force at time t = h as 

66 
»b=l 

-u(r ib )]. (28) 
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With Eq. tpS]) the velocities v*(/i) are obtained by 

Vi(h) = Vi(0) + — [F Cji (0) + F C)i (/i) 

+F f)i (0) + F f)i (/ i ) + 2F rii (/ i )] (29) 

Of course, at the same time one has to transfer also the force — Ff^(h) — F It i(h) to the 
appropriate nodes in the LB fluid. 

Our scheme is equivalent to the Heun algorithm 2 ^ which is a standard method for the 
solution of Langevin equation s 31 1 32 . If we set £o,i b — in Eqs. (J2~5|) and ()29j) the algorithm 
reduces to the velocity Verlet algorithm for the microcanonical ensemble 29 . 

V. RESULTS 

In this section we present several applications of our MD/LB method. First, we consider 
a neutral colloidal particle in a LB fluid to calculate its effective friction coefficient (part 
A) and to study long time tails in the translational and rotational velocity autocorrelation 
function (part B). In part C we consider a charged colloid in an electric field. In the following, 
we choose p = 1.0 mo/a 3 for the density and r] = up = 0.02381a 2 p/T for the shear viscosity 
of the LB fluid. 

A. Friction coefficient 

Consider a sphere with (hydrodynamic) radius that moves through a viscous fluid 
due to a gravitational field g. In the steady state, it experiences a drag force Fd which is 
proportional to its velocity U, according to Stokes law^ 3 -: 

F d = £U (30) 

where £ is the Stokes friction coefficient. In the case of no-slip boundary conditions, the 
friction coefficient is given by £ = QtitiR^ (with 77 being the shear viscosity of the fluid). In 
the steady state, the force F d is equal to the gravitational force on the sphere, 

F g = \kRI{p? ~ P)g, (31) 

where p p and p denote the density of the sphere and the fluid, respectively. Thus, from 
F g = Fd the friction coefficient £ can be determined by measuring the velocity U of the 
particle in the steady state. 
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As a convenient way to determine £ in a LB simulation, we follow here the scheme 
proposed by LadcU^. We consider a particle in a LB fluid. This system is put in a cubic box 
of volume V = L 3 with periodic boundary conditions in all three Cartesian directions x, y, z. 
The particle is held fixed by assigning an infinite mass to it. Then, a pressure gradient V x p 
is introduced in the x direction by applying a constant increment Aj x to the x component 
of the momentum density at each lattice node. In the steady state, the total force on the 
particle is balanced by the sum of the drag force F^ x = VAj x /r and the buoyancy force 
Fb = —~7rR 3 Aj x /T (remember that r denotes the time step of the LB simulation). Thus 
the Stokes friction coefficient is given by £l = F& iX /U x where the index L reflects that the 
particle is moving in a finite system of size L. Due to the periodic boundary conditions, £^ 
describes the friction of an array of spheres that sit on a cubic lattice with lattice constant L. 
An analytic expression for l/£z, in terms of an expansion of powers of 1/L was first derived 
by Hasimotc 34 assuming no-slip boundary conditions. It has the following form: 

1 Wl 2-837 + 4J. 1+ N 



£ L 67T77 \R h L L 3 
For our fluid-particle coupling scheme, we can only use this formula for high values of 
£o- Then, no-slip boundary conditions are approximately recovered as we shall see in the 
following. 

In Fig. El is plotted as a function of 1/L for different values of £ from £ = 3.3 m Q /r 
to £o — 19-8 m o/ T - The radius of the particle is R = 2.5a in this case. For 1/L < 0.04a _1 
the data can be well described by a linear 1/L dependence of the form 

f = r-4' (33) 

C.L t,oo L, 

where £oo denotes the friction coefficient for an unbounded system. Fits with Eq. f!33j) are 
shown in Fig. El as dashed lines. In these fits, the slope B changes only slightly with £o : We 
find B = 6.419ar/mo at £o — 3.3 ttiq/t and B = 6.438ar/mo at £ — 19-8 m o/ T - Moreover, 
for £ = 19.8 m /r the data over the full 1/L range can be well described by a fit with 
Eq. (}3*2|) (bold solid line in Fig. |2J). In the latter fit only the hydrodynamic radius R^ was 
used as a fit parameter for which we find R^ = 3.05a. The good quality of the fit indicates 
that at £o — 19.8 mo/r, no-slip boundary conditions are almost recovered whereas at smaller 
values of £o mixed stick-slip boundary conditions are obtained. This can be also inferred 
from the inset of Fig. El where we have plotted Rh/R as determined from the fitted ^ via 
R h /R = U/(Qn V R). 

12 



One may wonder why the hydro dynamic radius is about 20% higher than the assigned 
radius R = 2.5a of the particle. But this is just an artifact of the discrete nature of the LB 
fluid. This artifact can be reduced by increasing the size of the particle. 

The discrete nature of the LB fluid is also reflected in the dependence of £ on the size of 
the particle R. In Fig. |2l we show £ as a function of R for the two system sizes L = 60a and 
L = 80a. We see from the figure that, at small values of R, £l does not increase linearly 
with the sphere radius R and thus the ratio Rh/R is not a constant in this regime. Only, 
for R > 3.0a, £^ oc R seems to hold to a good approximation. 

We have seen that, with the frictional coupling scheme used in this work, no-slip boundary 
conditions at the surface of a colloidal particle can be nearly realized. For the system 
considered in this section, one has to choose a value of £o around 20 ttiq/t to obtain no-slip 
boundary conditions. Note that the limit of small values of £o is also of interest. As we shall 
demonstrate in a forthcoming publication^, the frictional coupling scheme can be also used 
to model walls at which mixed stick-slip boundary conditions hold. This is an important 
issue for the modeling of fluids in nanoscopic slits. 

B. Long time tails 

In this section, we consider again a neutral, spherical particle in a LB fluid. Our aim is 
to study the normalized translational and rotational velocity autocorrelation function which 
we denote by C v (t) and C u (t), respectively. These functions can be simply calculated by 
making use of linear response theory. To this end, we consider a colloidal particle in a LB 
fluid at rest. At t = 0, we give the particle a translational or rotational velocity. The decay 
of these velocities with time, normalized by the initial values, yields then the functions C v (t) 
and Cuj(t), respectively. 

Fig. 0] illustrates the short time decay of C v (t) for different values of £o- The initial 
decay of these functions is given by the exponential functions f(t) = exp (— ||t) which are 
shown as dotted lines in Fig. 0] (note that we have assigned the mass M = 120m to the 
particle). Thus, the short time behavior of C v (t) is completely controlled by the relaxation 
time tb = M/£o- 

As we can also see in Fig. for t > t-q the the average fluid velocity at the surface of 
the particle essentially equals the velocity of the particle. The dashed lines in Fig. 0] are 
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averaged velocities v s of the fluid at the boundary points of the particle (normalized by the 
initial velocity V(0) of the particle). Obviously, v s /V(0) matches the function C v (t) at a 
time which is of the order of tb- 

For long times one expects the occurrence of a long time tail in C v (t), i.e., a decay 
with the power law f(t) = At~ 3 / 2 — . The theoretical prediction for the prefactor is A = 
Mj (12p) [7ti/]~ 3 / 237 i 38 (with v = rj/p being the kinematic viscosity of the fluid). According 
to this prediction the prefactor A does not depend on the details of the particle such as its 
radius R. The physical origin of the long time tail is the conservation of momentum, which 
is transported away diffusively from the particle. Since the momentum transport in a fluid 
is spatially long-ranged, one expects the presence of finite size effects if the decay of velocity 
correlations of a colloidal particle in a finite simulation box is considered. 

In Fig. C v {t) is plotted for different system sizes L for £ = 6.6 m /r and, as a dotted 
line, the theoretically predicted long time tail is shown (note that no fit parameters are 
involved). We see that the theoretical result match perfectly with the simulation data. 
One can also infer from the figure the expected finite size effects. At a given L, C v seems 
to approach a constant at long enough times. This is just a consequence of momentum 
conservation: At long times the particle has completely transfered its initial velocity v to 
the LB fluid and then the whole system moves with a constant velocity Vo/N to t where iV to t 
is the total number of lattice nodes. 

Also the angular velocity correlation function C w (t) exhibits a long time tail but now the 
exponent for the power law decay is —5/2, i.e. f(t) = Bt~ 5 ^ 2 . The theoretical prediction for 
the prefactor is B = nl/p [Anu]^ 5 ^ 239 . In order to compare the latter theoretical prediction 
to the simulation result, one has to estimate the moment of inertia I of the rotating sphere. 
On the one hand, / has a contribution I\ from the shell of the particle that consists of 66 
boundary points. On the other hand, at long times the moment of inertia is also affected 
by the rotating fluid inside the sphere which leads to a second contribution I2 to I. 1\ can 
be estimated by I\ = ^MR 2 = 500 m a 2 (i.e., the value of a hollow sphere with radius R). 
We have computed I2 numerically by considering all the lattice nodes of mass mo inside the 
sphere to obtain the value I2 = 228 mod 2 . Thus, the total moment of inertia of the rotating 
sphere is / = Jx + 12 — 728 m a 2 . The theoretical prediction using this value of / is shown 
in Fig. |U1 in comparison with the simulation data for different values of L. Obviously, the 
theoretical prediction is in nice agreement with the numerical data. 
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As we have already mentioned, we use linear response theory to determine C v {t) and C w (£) 
by considering a kicked particle in a LB fluid which is initially at rest. Of course, we can also 
calculate the velocity correlation functions from thermal fluctuations of the translational and 
rotational velocity, respectively, and this should yield identical results. Indeed, this is the 
case as we can infer from Fig. [7|and Fig. |S] where C v (t) and C u (t), respectively, are shown 
for L = 40a and £o = 6.6 ttiq/t. 

We have seen in this section that we recover the the theoretical predictions for the long 
time tails in C v (t) and C w (t). This is in agreement with the study of Lobaskin et a/. 16 who 
use a slightly different particle-fluid coupling scheme (see above). 



C. Charged Colloids 

Now we demonstrate that our hybrid MD/LB scheme can be applied to charged colloidal 
systems. To this end, we consider a charged spherical particle (macroion) that is immersed 
in a fluid of small ions and a neutral "hydrodynamic background" which is modeled by a 
LB fluid. This system is then studied in an electric field to determine the drift velocity and 
the effective dynamic charge of the macroion. In the following, we first introduce the model 
and the simulation details, before we present the results of the simulation. 



1. Model and simulation details 



A potential that describes surprisingly well the effective interactions between macroions 
in a charged colloidal suspension is the Debye-Hiickel (DH) potential, which is the solution 
of the linearized Poisson-Boltzmann equation^. This potential has a Yukawa form, 

u{r) = if eXp( - r/AD) , (34) 
r 

where Ad is the so-called Debye screening length and K is a constant depending in partic- 
ular on the charge of the macroions Z m . Eq. (J53jl describes the screening of the Coulomb 
interaction between positively charged macroions due to the presence of small ions in the 
system, namely negatively charged counterions and positively charged coions. The screen- 
ing length is explicitly given by Ad = (47t/b zips) 1 / 2 where 1-q is the Bjerrum length and 
p s is the average density of microscopic ions of type s (either counterions or coions). In 
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principal the DH potential should only be valid for small charge and surface potential of 
the macroions. However, it turns out that also systems with highly charged particles can 
be rather well described by the DH potential if one introduces an effective charge for 
the macroions. One possibility to determine Z^ experimentally is via electrophoresis where 
the Z^ is extracted from the measurement of the electrophoretic mobility (see below). We 
shall address here the problem of estimating Z^ via electrophoresis by means of our LB/MD 
simulation method. 

To this end, we investigate a charged colloidal system in the framework of the primitive 
model 1 . Thus, we consider a system of a positively charged macroion of charge Z m and small 
ions of charge Z- x = — 1 (counterions) and of charge Z c = 1 (coions). Of course, the number 
of counterions and coions is chosen such that charge neutrality of the system holds. The 
interaction potential between a particle of type a and a particle of type (3 (a, (3 — m, i, c) 
separated by a distance r from each other is given by 

Z Z e 2 

u a /3 = -j— — _ h A af3 exp{-B af3 (r - a aP )} (35) 

47re r e r 

where e is the elementary charge and e r and eo are the reduced dielectric constant (which 
we set to e r = 80 for water at room temperature) and the vacuum dielectric constant, 
respectively. The parameter a a p is the distance between two ions at contact, a a p = R a + Rp, 
where R a is the radius of an ion of type a. In the following, we use R m = 20 A and R\ = R c = 
1 A. The exponential in Eq. ()35|) is an approximation to a hard sphere interaction for two ions 
at contact. For the parameters A a/3 we choose A mm = 1.84 eV, A mi = A mc = 0.0556544 eV, 
and An = A ic = A cc = 0.0051 eV. The parameters B a p are all set to 4.0 A -1 . 

We have done simulations for two different systems: The first system is a mixture of 
a macroion of charge Z m = 121 with 471 counterions and 350 coions. The second system 
contains a macroion of charge Z m = 255, 555 counterions, and 300 coions. In both cases, the 
ions are placed in a cubic simulation box of linear size L = 160 A using periodic boundary 
conditions. All the simulations are done at T = 297 K. The Debye screening length Ad is 
for both systems around 7.5 A. For the masses of the macroion and the small ions we have 
chosen 60 atomic units and 4 atomic units, respectively. Thus, the mass of the macroion 
is a factor 15 times the mass of a small ion. The Coulomb part of the potential and the 
forces was evaluated by means of Ewald sums with a constant a = 0.05 and by using for the 
Fourier part of the Ewald sum all k vectors of magnitude less than k c = 2n^/66/ L 29 ^ . 

16 



For the LB fluid to which the ions are coupled we use a cubic lattice with 40 3 lattice 
nodes. Hence, since the size of the simulation box is L = 160 A, the lattice constant of the 
LB fluid is a = 4.0 A. The counterions and coions are treated as point particles with respect 
to the LB fluid, i.e., each counterion and coion are equivalent to a sinlge boundary point on 
the surface of the colloid. The force on each of the ions is calculated by Eq. ((15)) . The value 
of the applied random force on the ion is calculated according to Eq. (}2*2*)) and is equal to the 
random force acting on a single boundary node. The macroion is seen by the LB fluid as a 
sphere of radius -R mj LB = 2.5a using the model with 66 boundary points and £ = 6.6m /r. 
The choice -R m ,LB < Rm is important if systems with more than one macroion are simulated. 
In this case the smaller -R mj LB prevents that two macroions get in close contact with respect 
to the LB fluid. 

The equations of motion were integrated with a time step of 1 fs. This very small time step 
is necessary because we consider explicitly counterions and coions as microscopic particles. 
A larger time step could be used if the exponential term in the potential, Eq. (JHHJ), is replaced 
by a softer repulsive potential. Moreover, for the update of the LB fields one could use a 
larger time step as for the MD part as it was done in Refs . iei24 . However, such optimizations 
were not necessary for the problems that are considered in the next section. 

2. A colloidal particle in an electric field 

The two systems with Z m = 121 and Z m = 255 are now studied in a constant electric 
field. This leads to a drift velocity of the macroion in the direction of the field (in the 
following we apply a field E x in x direction). In the linear response regime, Vd is linearly 
related to E x , Vd = nE x , where [i is the so-called electrophoretic mobility 1 . The latter 
quantity is of particular interest because, the experimental determination of fi allows to 
extract the effective charge of the macroio n 41 ' 42 ^ . In the following, we determine Vd as a 
function of E x and we estimate the effective charge of the macroion from the number of 
counterions that move with the macroion when the electric field is switched on. 

We first equilibrated our system for 25000 time steps without the electric field and without 
coupling the system to the LB fluid. Fig. EH is an equilibrium snapshot of the system with 
macroion charge Z m = 255. One can clearly identify a spherical region around the macroion 
(big sphere), the so-called Debye layer, that contains an excess of counterions (small light 
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gray spheres) whereas the coions (small black spheres) are almost excluded from this region. 
Far away from the Debye layer the small ions are randomly distributed. The presence of the 
Debye layer can be also inferred from the radial density distributions p c of counterions and 
coions around the macroion which are shown in Fig. El for the two systems with different 
macroion charges. In the inset of Fig. El the instantaneous temperature of the small ions 
is shown during a simulation of 10000 time step. We see that it fluctuates correctly around 
the assigned value for the temperature, T = 297 K. 

After the equilibration, the system was coupled to the LB fluid and the electric field E x 
was switched on. We did runs for several values of E x ranging from E x = 0.0025 V/A to 
E x = 0.02 V/A. For each value of the electric field, runs over 300000 to 500000 steps were 
done. After reaching the steady state within about 10000 time steps, positions and velocities 
of the ions were stored every 500 steps to obtain the averaged quantities that are presented 
in the following. 

Due to the electric field E x the macroion and the coions move in the positive x direction 
whereas the counterions experience a force in the opposite direction. This leads to a dis- 
tortion of the Debye layer which gives rise to a force opposing the motion of the macroion. 
This retarding force coupled with the viscous drag due to the fluid balances the force on the 
macroion due to E x in the steady state. Fig. ITT1 shows steady state configurations of the ions 
for the system with macroion charge Z m = 255 for two different choices of E x , namely for 
E x = 0.01 V/A (top) and E x = 0.02 V/A (bottom), where the positive x direction is marked 
by a black arrow. It can be seen from the picture that there are more counterions behind the 
macroion than in front of it leading to the distorted counterion charge distribution. A high 
value of the electric field (Fig. ^2 bottom) strips off the counterions from the shell moving 
along with the macroion and hence the counterion cloud becomes more diffuse. 

In order to obtain an estimate of the effective charge of the macroion, we determine the 
average number of counterions N* that move along with the macroion, i.e. those counterions 
are counted that have a positive drift velocity in x direction. In the limit of small electric field 
E x , i.e. in the linear response regime, \Zi\N* (here Zi = —1) is a measure of the effective 
charge Z^ of the macroion. Hence, the ratio N*/Z m approaches one towards vanishing 
values of E x if Z* m is equal to Z m . Fig. fT2l shows N*/Z m as a function of E x for the two 
different systems with Z m = 121 and Z m = 255. In both cases, the linear response regime is 
obviously not approached, even for small considered values of E x . However, an extrapolation 



18 



of the two curves in Fig. ^] to zero electric fields yields values close to one and thus, we 
can conclude from the data that is close to the bare charge Z m for the systems under 
consideration. We also see in Fig. that the slope for small E x is smaller for the system 
with Z m = 255 than for the one with Z m = 121. This finding shows that, the smaller Z m 
is, the smaller electric fields are required to obtain a reasonable estimate of the effective 
charge Z^. That the linear response regime is not reached for our smallest considered values 
of E x , is illustrated in the inset of Fig. ^] which shows the drift velocity Vd as a function 
of E x . The reason that we did not do simulations for lower values of E x is that the drift 
velocity approaches the order of the thermal velocity then and thus, it is difficult to yield a 
reasonable statistics. 

The finding Z^ « Z m is in agreement with recent MD simulation o 43 ' 44 ' 45 of similar systems 
where the effective charge Z^ was estimated from the potential of mean force between the 
macroions. 



VI. CONCLUSIONS 

We have presented a hybrid MD /LB method for the simulation of colloidal systems that 
has been applied to several simple colloidal systems. Our method can be seen as an alter- 
native to Ladd's coupling schem o 12 ! 13 as well as to the one proposed by Lobaskin et a/ . 16 i 24 
We think that all these methods (including ours) have their advantages and disadvantages 
and it depends on the problem whether one might prefer one or the other method. 

We have applied our method to the simulation of charged colloidal systems where, apart 
from macroions, counterions and coions are considered in the framework of the primitive 
model. For such systems, our simulation technique has been used to get insight into the 
properties of colloids in an external electric field. This is of particular interest for the 
understanding of experimental studies where one extracts effective charges of macroions 
from the measured electrophoretic mobilit y 41142143 . A more extensive study on this issue is 
in preparation. 

The model with the explicit consideration of small ions has also a drawback which was 
already pointed out in Ref.— : Only relatively small size ratios between macroions and small 
ions are accessible. Furthermore, only relatively small systems can be simulated due to the 
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long-ranged character of the electrostatic interactions. Thus, one looses partly the advan- 
tage of the LB method that mesoscopic length and time scales can be covered. However, 
as we have mentioned above, the DH potential describes the effective interactions between 
macroions surprisingly well. Therefore, we plan to study systems of, say, 1000 macroions 
which interact with each other via an effective DH potential. This allows then the investi- 
gation how the long-time diffusion of systems of charged colloidal particles at intermediate 
densities is affected by hydrodynamic interactions. 
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VII. LIST OF CAPTIONS 



Fig. ^ Sketch of the model for a colloidal particle. The centres of the small spheres 
represent the points at which the colloidal particle interacts with the LB fluid. The cylinders 
that connect the spheres with each other are just guides to the eye. 

Fig. El for particle of radius R = 2.5a as a function of 1/L for the indicated values 
of £o- The dashed lines are fits with Eq. ()33|) and the solid lines is a fit with Eq. (}3"2*|) for 
£o = 19.8mo/V using the hydrodynamic radius i?h as a fit parameter. The inset shows the 
variation of the hydrodynamic radius normalized by the assigned sphere radius R as a 
function of £ (see text). 

Fig- El £l as a function of the particle radius R for £ = 13.2m /r and the system sizes 
L = 60a and L = 80a. 

Fig. HI Time-dependent velocity correlation function C v (t) (solid lines) of the colloidal 
particle and normalized average velocity v s (t)/V(0) of the LB fluid at the surface of the 
particle (dashed lines) for a) £o = 0.66 mo/r, b) £o = 3.3 m /r, c) £ = 6.6 m /r, and d) 
£o = 13.2 m /r. The box length was set to L — 40a, the radius of the particle to R = 2.5a, 
and its mass to M = 120m . The dotted lines are exponential functions f(t) = exp (— f|i) 
that describe the short-time decay of C v (t) (note that no fit parameter is involved). 

Fig. El C v (t) for £ = 6.6 m /r at the indicated system sizes. The dotted line shows the 
power law f(t) = At~ 3 ^ 2 (see text). 

Fig. El The same as in Fig. El but now for the angular velocity correlation function C u {t) 
and for £ = 13.2mo/T. The dotted line shows the power law fit) = Bt~ 5 ^ 2 (see text). 

Fig. El Velocity autocorrelation function for £o = 6.6mo/r and L = 40a as calculated for 
a kicked particle (solid line) and with thermal fluctuations (circles) . 

Fig. El The same as in Fig. [7| but now for the angular velocity correlation function. 

Fig. El Representative configuration of counterion and coion distribution around a 
macroion of charge Z m = 255. The macroion is the central large grey sphere, the 555 
counterions are the light gray coloured small spheres and the 300 coions are the black small 
spheres. 

Fig. m3 Plot of the radial density distribution p m (r) of counterions and coions around 
the macroion for the two indicated charges Z m . The system with Z m = 255 contains 555 
counterions and 300 coions whereas there are 471 counterions and 350 coions in the system 
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with Z m = 121. The inset shows the fluctuation of the temperature T of the counterions 
around its assigned value (solid line). 

Fig. ^2 Snapshot of the ionic distribution around the macroion (for the system with Z m = 
255) with an applied electric field of E x = 0.01 V/A (top) and E x = 0.02 V/A (bottom), 
respectively. The electric field is applied in the direction of the black arrow. The rest is 
similar to Fig. El 

Fig. El N*/Z m as a function of the external electric field E x . N* is the average number of 
counterions that move with the particle. The inset shows the drift velocity v& as a function 
of the electric field E x . 
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